
/*************
author:: libo
license:: MIT
*****************/
#include<Mesh/_Func_.h>
#include<Mesh/libcell_iterator.h>
#include<tools_vector.h>
#define iterator_v lib_cell_iterator_vertices
#define iterator_f lib_cell_iterator_faces
#define iterator_c lib_cell_iterator_cells
#define iterator_edge lib_cell_iterator_edge

#define quote lib_cell_quote

LB_Matrix* mesh_compute_antisymmetric_tensormatrix(double**p,LB_Matrix* mat)
{
    double *data=(double*)(mat->data);
    int rows=mat->rows,cols=mat->cols;
    if(rows>cols){return NULL;}
    for(int i=0;i<rows;i++)
    {
        for(int j=0;j<cols;j++)
        {
            data[i*cols+j]=p[i+1][j]-p[0][j];
        }
    }
    return mat;
}

static template_v* Mesh_get_vertexp(struct Mesh *m,int i)
{
    return (template_v*)(m->vertices->find(m->vertices,i));
}
static template_c*  Mesh_get_cellp(struct Mesh*m,int i)
{
    return (template_c*)(m->cells->find(m->cells,i));
}
static template_f* Mesh_get_facep(struct Mesh*m,int i)
{
    return (template_f*)(m->faces->find(m->faces,i));
}
static template_hf* Mesh_get_halffacep(struct Mesh*m,int i)
{
    return (template_hf*)(m->halffaces->find(m->halffaces,i));
}
static template_v Mesh_get_vertex(struct Mesh* own,int id)
{
    template_v *v=Mesh_get_vertexp(own,id);
    if(v!=NULL)
    {
        return *v;
    }
    template_v v1;
    Vertex_init_(&v1);
    return v1;
}
static template_c Mesh_get_cell(struct Mesh* own,int id)
{
    template_c* c=Mesh_get_cellp(own,id);
    if(c!=NULL)
    {
        return *c;
    }
    template_c c1;
    Cell_init_(&c1);
    return c1;
}

static template_f Mesh_get_face(struct Mesh*own,int id)
{
    template_f *f=Mesh_get_facep(own,id);
    if(f!=NULL)
    {
        return *f;
    }
    template_f f1;
    Face_init_(&f1);
    return f1;
}

static template_hf Mesh_get_halfface(struct Mesh*own,int id)
{
    template_hf* hf=Mesh_get_halffacep(own,id);
    if(hf!=NULL)
    {
        return *hf;
    }
    template_hf hf1;
    HalfFace_init_(&hf1);
    return hf1;

}
template_f* Mesh_get_facev(struct Mesh*own,template_v** temp_v,int size)
{
    // Node* node=temp_v[0]->faces;
    // bool b=false;
    // while(node!=NULL)
    // {   b=true;
    //     for(int i=1;i<size;i++)
    //     {
    //         if(temp_v[i]->faces==NULL){return NULL;}
    //         if(node_find(temp_v[i]->faces,node->value)==NULL)
    //         {
    //             b=false;
    //             break;
    //         }
    //     }
    //     if(b)
    //     {
    //         return (template_f*)(node->value);
    //     }
    //     node=(Node*)(node->Next);
    // }
    // return NULL;
    bool b=false;
    for(int i=0;i<temp_v[0]->facesn.size;i++)
    {
        b=true;
        template_f* ref=(template_f*)(temp_v[0]->facesn.extra_data[i]);
        for(int j=1;j<size;j++)
        {
            if(lb_vector1_find_voidptr(&(temp_v[j]->facesn),ref)==-1)
            {
                b=false;
                break;
            }
        }
        if(b)
        {
            return ref;
        }
    }
    return NULL;
}

static template_c* Mesh_get_cellv(struct Mesh*own,template_v** temp_v,int size)
{
    // Node* node=temp_v[0]->cells;
    // bool b=false;
    // while(node!=NULL)
    // {
    //     b=true;
    //     for(int i=1;i<size;i++)
    //     {
    //         if(temp_v[i]->cells==NULL){return NULL;}
    //         if(node_find(temp_v[i]->cells,node->value)==NULL)
    //         {
    //             b=false;
    //             break;
    //         }
    //     }
    //     if(b)
    //     {
    //         return (template_c*)(node->value);
    //     }
    //     node=(Node*)(node->Next);
    // }

    // return NULL;

    bool b=false;
    for(int i=0;i<temp_v[0]->cellsn.size;i++)
    {
        b=true;
        template_c* rec=(template_c*)(temp_v[0]->cellsn.extra_data[i]);
        for(int j=1;j<size;j++)
        {
            if(lb_vector1_find_voidptr(&(temp_v[j]->cellsn),rec)==-1)
            {
                b=false;
                break;
            }
        }
        if(b)
        {
            return rec;
        }
    }
    return NULL;

}
static template_c* Mesh_get_cellf(struct Mesh* own,template_f** temp_f,int size)
{
    bool b=false;
    int i=0;
    for(i=0;i<2;i++)
    {
        if(temp_f[0]->halffaces[i]->cell==NULL)
        {
            continue;
        }
        b=true;
        for(int j=1;j<size;j++)
        {
            if(temp_f[j]->halffaces[0]->cell==temp_f[0]->halffaces[i]->cell||temp_f[j]->halffaces[1]->cell==temp_f[0]->halffaces[i]->cell)
            {
                continue;
            }
            else
            {
                b=false;
                break;
            }
        }
        if(b)
        {
            return temp_f[0]->halffaces[i]->cell;
        }
    }
    return NULL;
}

static int Mesh_num_c(struct Mesh* own)
{
    return own->cells->size;
    //return own->cells.size();
}
static int Mesh_num_f(struct Mesh* own)
{
    return own->faces->size;
    //return own->faces.size();
}
static int Mesh_num_v(struct Mesh* own)
{
    return own->vertices->size;
	//return own->vertices.size();
}
static int Mesh_num_hf(struct Mesh* own)
{

    return own->halffaces->size;
    //return own->halffaces.size();
}
static template_v* Mesh_create_vertex(struct Mesh* own)
{
    template_v* v=(template_v*)malloc(sizeof(template_v));

    Vertex_init_(v);
    v->id=own->vertex_id++;
    //own->vertices[v->id]=v;
    own->vertices->insert(own->vertices,v->id,v);
    return v;
}
static template_v* Mesh_create_vertexv(struct Mesh*own,double* p,int size)
{

    template_v*v= Mesh_create_vertex(own);
    v->point_size=size;
    v->point=(double*)malloc(sizeof(double)*size);
    for(int i=0;i<size;i++)
    {
        v->point[i]=p[i];

    }
    return v;
}
static template_f * Mesh_create_face(struct Mesh* own)
{
    template_f* f=(template_f*)malloc(sizeof(template_f));
    Face_init_(f);
    f->halffaces[0]->cell=NULL;f->halffaces[1]->cell=NULL;
    f->id=own->face_id++;

    //own->faces[f->id]=f;
    own->faces->insert(own->faces,f->id,f);
    return f;

}
static template_c* Mesh_create_cell(struct Mesh* own)
{
    template_c *c=(template_c*)malloc(sizeof(template_c));

    Cell_init_(c);
    c->id=own->cell_id++;

    //own->cells[c->id]=c;
    own->cells->insert(own->cells,c->id,c);

    return c;
}
//如果返回NULL,意味创建失败
template_hf * Mesh_create_halfface(struct Mesh* own,template_f* f,template_v** temp_v,int size)
{
    if(f==NULL){return 0;}
    if(f->halffaces[0]->cell!=NULL&&f->halffaces[1]->cell!=NULL)
    {
        return NULL;
    }


    template_hf *hf=NULL;
    if(f->halffaces[0]->vertices==NULL)
    {
        hf=f->halffaces[0];
    }
    else if(f->halffaces[1]->vertices==NULL)
    {

        hf=f->halffaces[1];
    }
    else
    {
        if(f->halffaces[0]->cell==NULL)
        {
            hf=f->halffaces[0];
        }
        else if(f->halffaces[1]->cell==NULL)
        {
            hf=f->halffaces[1];
        }
        else
        {
            printf("流形的网格拓扑错误\n");
            for(int i=0;i<size;i++)
            {
                printf("%d ",temp_v[i]->id);
            }
            printf("\n");
            return NULL;
        }
    }
    hf->id=own->halfface_id++;


//#ifndef LIBCELL_MEMORY_SAVING_MODE
    own->halffaces->insert(own->halffaces,hf->id,hf);

    //SAFE_FREE(hf->vertices);
    if(hf->vertices==NULL)
    {
        hf->vertices=(template_v**)malloc(sizeof(template_v*)*size);
    }
    memmove(hf->vertices,temp_v,sizeof(template_v*)*size);
    // for(int i=0;i<size;i++)
    // {
    //     hf->vertices[i]=temp_v[i];
    // }
    hf->vertices_size=size;


    return hf;
}
//以下程序表明点到face的遍历没有不能保证顺(逆)时针(当cell是曲面)
template_f* Mesh_create_facev(struct Mesh*own,template_v** temp_v,int size)
{
	template_f*f=Mesh_get_facev(own,temp_v,size);
	if(f==NULL)
	{
	    f=Mesh_create_face(own);
	}
	else
	{
	    return f;
	}
	f->vertices_size=size;
	f->vertices=(template_v**)malloc(sizeof(template_v*)*size);
    for(int i=0;i<size;i++)
    {

        f->vertices[i]=temp_v[i];
        // temp_v[i]->faces=node_overlying(temp_v[i]->faces,(void*)f);
        lb_vector1_pushback_voidptr(&(temp_v[i]->facesn),f);

    }
    return f;
}
//要判断是否已经创建相应的cell
template_c* Mesh_create_cellf(struct Mesh* own,template_hf**temp_hf,int size)
{
//#ifdef __DEBUG__
    for(int i=0;i<size;i++)
    {
        if(temp_hf[i]==NULL)
        {
            if(own->manifold_require==1)
            {
                for(int j=0;j<size;j++ )
                {
                    if(temp_hf[j]!=NULL)
                    {
                        template_f* f=temp_hf[j]->face;
                        //
                        Mesh_delete_halfface(own,*(temp_hf[j]),true);
                        template_hf*hf=(f->halffaces[0]==temp_hf[j]?f->halffaces[1]:f->halffaces[0]);
                        if(Mesh_get_halffacep(own,hf->id)==NULL&&f->halffaces[0]->cell==NULL&&
                            f->halffaces[1]->cell==NULL)
                        {
                            for(int k=0;k<f->vertices_size;k++)
                            {
                                Vertex_remove_f(f->vertices[k],f);
                            }
                            Mesh_delete_face(own,*f,true);
                            free_Face(f);
                        }
                    }
                }
            }
            return NULL;
        }
    }
//#endif

    if(temp_hf[0]->cell!=NULL)
    {
        //exist cell
        return temp_hf[0]->cell;
    }
    template_c* c=Mesh_create_cell(own);
    if(c==NULL){printf("c is null\n");return NULL;}

    // Node*node_v=NULL;
    LB_Vector_voidptr vec;
    lb_vector_init_voidptr(&vec);
    c->halffaces_size=size;
    c->halffaces=(template_hf**)malloc(sizeof(template_hf*)*size);
    for(int i=0;i<size;i++)
    {
        temp_hf[i]->cell=c;
        c->halffaces[i]=temp_hf[i];
        //node=node_pushback(node,(void*)temp_hf[i]);
        // printf("%d ",temp_hf[i]->vertices_size);
        for(int j=0;j<temp_hf[i]->vertices_size;j++)
        {
            // if(node_findr(node_v,temp_hf[i]->vertices[j])==NULL)
            if(vec.find(&vec,temp_hf[i]->vertices[j])==-1)
            {
                vec.push_back(&vec,temp_hf[i]->vertices[j]);
                // printf("pushback:%d  ",temp_hf[i]->vertices[j]->id);
                // node_v=node_pushback(node_v,temp_hf[i]->vertices[j]);
            }
        }
        // printf("\n");
    }

    //c->halffaces=node_reverse(node);

    // node_v=node_reverse(node_v);
    //Node* temp_nodev=node_v;
    // printf("end mesh create cell1%d\n",node_size(node_v));
    // int le=node_size(node_v);
    int le=vec.size;
    // printf("le :%d\n",le);
    c->vertices=(template_v**)malloc(sizeof(template_v*)*le);
    c->vertices_size=le;
    // int i=0;
    // for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    for(int i=0;i<vec.size;i++)
    {
        c->vertices[i]=(template_v*)(*(vec.at(&vec,i)));
        lb_vector1_pushback_voidptr(&(c->vertices[i]->cellsn),c);
        // c->vertices[i]=(template_v*)(nit->value);
        // c->vertices[i]->cells=node_overlying(c->vertices[i]->cells,(void*)c);
        // i++;
    }
    vec.clear(&vec);
    // free_node(node_v);
    return c;
}
//应该写 Mesh_get_cellv，并用在下面判断(但速度会变慢)
//以下程序也不能保证点到cell的遍历是顺(逆)时针(当cell表示面时)
//创建单形使用的函数
//如果返回NULL意味创建失败
template_c* Mesh_create_cellv(struct Mesh* own,template_v** temp_v,int size)
{
    template_c* ccc= Mesh_get_cellv(own,temp_v,size);

    if(ccc!=NULL)
    {
        return ccc;
    }

    template_hf** hfs=(template_hf**)malloc(sizeof(template_hf*)*size);
    LB_Vector_voidptr vec;
    lb_vector_init_voidptr(&vec);

    // template_v** temp_v1=(template_v**)malloc(sizeof(template_v*)*(size-1));

    for(int i=0;i<size;i++)
    {
        int temp_i=0;
        vec.size=0;
        for(int j=0;j<size;j++)
        {
            if(j!=i)
            {
                vec.push_back(&vec,temp_v[j]);
                // temp_v1[temp_i]=temp_v[j];
                // temp_i++;
            }
        }
        template_v**temp_v1=(template_v**)(vec.extra_data==NULL?vec.data:vec.extra_data);
        if((size-i-1)%2==0)
        {}
        else if(size>2)
        {
            template_v* vv=temp_v1[0];
            temp_v1[0]=temp_v1[1];
            temp_v1[1]=vv;
        }

        template_f* f=Mesh_create_facev(own,temp_v1,size-1);
        hfs[i]=Mesh_create_halfface(own,f,temp_v1,size-1);

        if(hfs[i]==NULL)
        {
            printf("create cellv is NULL\n");
            for(int j=0;j<size-1;j++)
            {
                printf("v : %d \n",temp_v1[j]->id );
            }
            printf("f :%d\n",f->id);
            if(own->manifold_require==1)
            {
                Int_RB_Tree* temp_tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
                int_rb_tree_init(temp_tree);
                for(int j=0;j<i;j++)
                {
                    own->halffaces->erase(own->halffaces,hfs[j]->id);
                    SAFE_FREE(hfs[j]->vertices);
                    hfs[j]->vertices_size=0;
                    if(hfs[j]->face->halffaces[0]->cell==NULL&&hfs[j]->face->halffaces[1]->cell==NULL)
                    {
                        temp_tree->insert(temp_tree,hfs[j]->face->id,hfs[j]->face);
                    }
                }
                if(f->halffaces[0]->cell==NULL&&f->halffaces[1]->cell==NULL)
                {
                    temp_tree->insert(temp_tree,f->id,f);
                }
                for(auto it=temp_tree->begin(temp_tree);it.it!=NULL;it++)
                {
                    template_f* ff1=(template_f*)(it.second);
                    own->faces->erase(own->faces,ff1->id);
                    for(int k=0;k<ff1->vertices_size;k++)
                    {
                        Vertex_remove_f(ff1->vertices[k],ff1);
                    }
                    free_Face(ff1);
                }
                int_rb_tree_free(temp_tree);
            }
            //free(temp_v1);
            free(hfs); vec.clear(&vec);
            return NULL;
        }
    }
    //free(temp_v1);
    vec.clear(&vec);
    template_c* c=Mesh_create_cell(own);
    c->vertices=(template_v**)malloc(sizeof(template_v*)*size);
    c->vertices_size=size;
    c->halffaces_size=size;
    c->halffaces=hfs;
    for(int i=0;i<size;i++)
    {
        c->vertices[i]=temp_v[i];

        lb_vector1_pushback_voidptr(&(temp_v[i]->cellsn),c);
        // temp_v[i]->cells=node_overlying(temp_v[i]->cells,(void*)c);
//#ifdef SIMPLEX_REQUIRE

        hfs[i]->cell=c;
        //c->halffaces=node_overlying(c->halffaces,(void*)(hfs[i]));
//#endif
    }
    //free(hfs);


    return c;
}

//
int Mesh_delete_vertex(struct Mesh*own,template_v v_,int burning)
{
    template_v*v=(template_v*)(own->vertices->find(own->vertices,v_.id));

    if(v==NULL)
    {
        return 0;
    }
    own->vertices->erase(own->vertices,v_.id);
    template_c*c;

    LB_Vector_voidptr vec;
    lb_vector_init_voidptr(&vec);
    // for(Node* nit=v->cells;nit!=NULL;nit=(Node*)(nit->Next))
    for(int i=0;i<v->cellsn.size;i++)
    {
        vec.push_back(&vec,v->cellsn.extra_data[i]);
    }

    for(int i=0;i<vec.size;i++)
    {
        template_c*c=(template_c*)(*(vec.at(&vec,i)));
        Mesh_delete_cell(own,*c,true);
    }
    vec.clear(&vec);

	// Node* node=node_copy(v->cells),*node1=NULL,*node2=NULL;

	// node1=node;
	// while(node1!=NULL)
	// {
 //        node2=(Node*)(node1->Next);
	//     c=(template_c*)(node1->value);
	//     Mesh_delete_cell(own,*c,true);
	//     node1=node2;

 //    }
	// free_node(node);


//#ifndef MANIFOLD_REQUIRE
    if(own->manifold_require==0)
    {
        // for(Node* nit=v->faces;nit!=NULL;nit=(Node*)(nit->Next))
        for(int i=0;i<v->facesn.size;i++)
        {
            vec.push_back(&vec,v->facesn.extra_data[i] );
            // vec.push_back(&vec,nit->value);
        }
        for(int i=0;i<vec.size;i++)
        {
            template_f*f=(template_f*)(*(vec.at(&vec,i)));
            Mesh_delete_face(own,*f,true);
        }
        vec.clear(&vec);
        // node=node_copy(v->faces);
        // template_f*f;
        // node1=node;
        // while(node1!=NULL)
        // {
        //     node2=(Node*)(node1->Next);
        //     f=(template_f*)(node1->value);
        //     Mesh_delete_face(own,*f,true);
        //     node1=node2;
        // }
        // free_node(node);
    }
//#endif

    if(burning)
    {
 	    free_Vertex(v);
	}

    return 1;
}
int Mesh_delete_halfface(struct Mesh*own,template_hf f_,int burning)
{
    template_hf*hf=(template_hf*)(own->halffaces->find(own->halffaces,f_.id));
    if(hf==NULL)
    {
        return 0;
    }

    if(hf->cell!=NULL)
    {
        Mesh_delete_cell(own,*((template_c*)(hf->cell)),burning);
        hf->cell=NULL;
    }
//#ifndef MANIFOLD_REQUIRE
    if(own->manifold_require==0)
    {
        if(hf->vertices!=NULL)
        {
            free(hf->vertices);
            hf->vertices=NULL;
        }
        hf->vertices_size=0;
        //own->halffaces.erase(it);
        own->halffaces->erase(own->halffaces,f_.id);
        //int_rb_tree_delete_node(own->halffaces,f_.id);
    }

//#endif

    return 1;
}
int Mesh_delete_face(struct Mesh*own,template_f f_,int burning)
{
    template_f*f=(template_f*)(own->faces->find(own->faces,f_.id));
    if(f==NULL)
    {
        return 0;
    }

    for(int i=0;i<2;i++)
    {
        if(f->halffaces[i]->cell!=NULL)
        {
            Mesh_delete_cell(own,*(f->halffaces[i]->cell),burning);
        }
    }

//#ifndef MANIFOLD_REQUIRE
    if(own->manifold_require==0)
    {
        own->faces->erase(own->faces,f_.id);

        //int_rb_tree_delete_node(own->faces,it);

        for(int i=0;i<f->vertices_size;i++)
        {
            Vertex_remove_f((template_v*)f->vertices[i],f);
        }
        if(burning)
        {
            free_Face(f);
        }
    }

//#endif
    return 1;
}
int Mesh_delete_cell(struct Mesh* own,template_c c_,int burning)
{
    template_c* c=(template_c*)(own->cells->find(own->cells,c_.id));
    if(c==NULL){return 0;}
    own->cells->erase(own->cells,c_.id);


    //handle with face
    template_f*f;
    for(int i=0;i<c->halffaces_size;i++)
    {
        Halfface_remove_c(c->halffaces[i],c);
        own->halffaces->erase(own->halffaces,c->halffaces[i]->id);
    }
    // Node* node=c->halffaces,*node1=NULL;
    // while(node!=NULL)
    // {
    //     node1=(Node*)(node->Next);
    //   //  f=(template_f*)(((template_hf*)node->value)->face);
    //     //face要甩掉cell
    //     //Face_remove_c(f,c);
    //     Halfface_remove_c((template_hf*)(node->value),c);

    //     // auto hfit=own->halffaces.find(((template_hf*)(node->value))->id);
    //     // if(hfit!=own->halffaces.end())
    //     // {
    //     //     own->halffaces.erase(hfit);
    //     // }
    //     own->halffaces->erase(own->halffaces,((template_hf*)(node->value))->id);
    //     node=node1;
    // }
    template_v*v;
    for(int i=0;i<c->vertices_size;i++)
    {
        v=(template_v*)c->vertices[i];
        Vertex_remove_c(v,c);
    }
//#ifdef MANIFOLD_REQUIRE
    if(own->manifold_require==1)
    {
        //Node* node_f=c->halffaces;
        //while(node_f!=NULL)
        for(int i=0;i<c->halffaces_size;i++)
        {
            f=c->halffaces[i]->face;
            if(f->halffaces[0]->cell==NULL&&f->halffaces[1]->cell==NULL)
            {
                for(int i=0;i<f->vertices_size;i++)
                {
                    Vertex_remove_f((template_v*)f->vertices[i],f);

                }

                own->halffaces->erase(own->halffaces,f->halffaces[0]->id);
                own->halffaces->erase(own->halffaces,f->halffaces[1]->id);
                own->faces->erase(own->faces,f->id);
                if(burning)
                {
                    free_Face(f);
                }
            }
            //node_f=(Node*)(node_f->Next);

        }
    //
    //如果要严格删除孤立点，则取消注释一下代码
        //*************
/*    for(int i=0;i<c->vertices_size;i++)
    {
        v=(template_v*)(c->vertices[i]);
        if(v->cells==NULL)
        {
            auto vit=own->vertices.find(v->id);
            if(vit!=own->vertices.end())
            {
                own->vertices.erase(vit);
            }
            free_Vertex(v);
        }

    }
  */
    }

//#endif

    if(burning)
    {
        free_Cell(c);

    }
    return 1;
}

static int Mesh_nface_is_boundary(struct Mesh* own,template_f f)
{
    if(f.halffaces[0]->cell==NULL||f.halffaces[1]->cell==NULL)
    {   return true;}
    return false;
}
static int Mesh_face_is_boundary(struct Mesh* own,template_f *f)
{

    return Mesh_nface_is_boundary(own,*f);
}

static int Mesh_cell_is_boundary(struct Mesh*own,template_c c )
{
    for(auto chfit=own->chf_begin(own,c);chfit!=own->chf_end(own,c);chfit++)
    {
        template_hf*hf=quote(chfit);
        if(Mesh_nface_is_boundary(own,*(hf->face)))
        {
            return 1;
        }
    }
    return 0;
}

//孤立点也是边界点
static int Mesh_vertex_is_boundary(struct Mesh* own,template_v v)
{
    if(v.facesn.size<=0)
    {
        return true;
    }
    // for(int i=0;i<v.facesn.size;i++)
    for(auto vfit=own->vf_begin(own,v);vfit!=own->vf_end(own,v);vfit++)
    {
        // if(Mesh_nface_is_boundary(own,*((template_f*)(v.facesn.extra_data[i]))))
        if(Mesh_nface_is_boundary(own,*vfit))
        {   return true;}
    }
    return false;
    // Node* node=v.faces;
    // if(node==NULL)
    // {
    //     return true;
    // }
    // while(node!=NULL)
    // {
    //     if(Mesh_nface_is_boundary(own,*((template_f*)(node->value))))
    //     {   return true;}
    //     node=(Node*)(node->Next);
    // }
    // return false;
}
int Mesh_face_is_in_cell(template_f*f ,template_c* c)
{
    for(int i=0;i<c->halffaces_size;i++)
    {
        //template_hf* hf=(template_hf*)(nit->value);
        if(c->halffaces[i]->face==f)
        {
            return true;
        }
    }
    return false;
}
int Mesh_vertex_is_in_cell(template_v*v,template_c*c)
{
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]==v)
        {
            return true;
        }
    }
    return false;
}
int Mesh_vertex_is_in_face(template_v*v,template_f*f)
{
    for(int i=0;i<f->vertices_size;i++)
    {
        if(f->vertices[i]==v)
        {
            return true;
        }
    }
    return false;
}
//
//给它一个初始面，他会返回boundary的node链.这种寻找边界的方式.
//返回半边
//用红黑树RB_Tree做标记，代替prop
//不确保保持流形

static Node* Mesh_node_of_boundary_face(struct Mesh* own,template_f *f_)
{
    Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree);

    if(f_==NULL)
    {
        int_rb_tree_free(tree);
        return NULL;
    }
    Node*node3=NULL;
    // Node* node=NULL;
    LB_Vector_voidptr vec;
    lb_vector_init_voidptr(&vec);
    template_f *f=NULL;

    template_hf* hf=NULL;
    f=(template_f*)(own->faces->find(own->faces,f_->id));
    if(f==NULL)
    {
        return NULL;
    }

    if(!Mesh_nface_is_boundary(own,*f_))
    {return NULL;}
    else
    {
        // node=node_overlying(node,(void*)f);
        vec.push_back(&vec,f);
        tree->insert(tree,f->id,f);
        //printf(" %d %d\n",f->vertices[0]->id,f->vertices[1]->id);
        if(f->halffaces[0]->cell==NULL)
        {
            hf=f->halffaces[0];
            //
            node3=node_overlying(node3,(void*)hf);
        }
        else if(f->halffaces[1]->cell==NULL)
        {
            hf=f->halffaces[1];
            node3=node_overlying(node3,(void*)hf);
        }

    }

   // 深度优先的遍历
//   printf("shendu\n");

    template_v *v=NULL;
    // Node* node2=NULL;
    //
    while(vec.size>0)
    // while(node!=NULL)
    {
        push:
        // f=(template_f*)(node->value);
        f=(template_f*)((vec.at(&vec,(vec.size-1)))[0]);
//        printf("f : %d\n",f->id);
        f->halffaces[0]->cell!=NULL?(hf=f->halffaces[0]):(hf=f->halffaces[1]);
        //tree->insert(tree,f->id,f);
        //node=node_remove(node,node);
        for(int i=hf->vertices_size-1;i>=0;i--)
        {
            v=(template_v*)(hf->vertices[i]);
            // for(Node* node2=v->faces;node2!=NULL;node2=(Node*)(node2->Next))
            for(int j=0;j<v->facesn.size;j++)
            {
                f=(template_f*)(v->facesn.extra_data[j]);
                // f=(template_f*)(node2->value);
                if(Mesh_nface_is_boundary(own,*f))
                {
                    void * var=tree->find(tree,f->id);
                    //if(rbt1==NULL)
                    if(var==NULL)
                    {
                        // node=node_overlying(node,f);
                        //printf(" %d %d\n",f->vertices[0]->id,f->vertices[1]->id);
                        vec.push_back(&vec,f);
                        void* temp_value=NULL;
                        if(f->halffaces[0]->cell==NULL)
                        {
                            temp_value=(void*)(f->halffaces[0]);
                            node3=node_overlying(node3,temp_value);
                        }
                        else if(f->halffaces[1]->cell==NULL)
                        {
                            temp_value=(void*)(f->halffaces[1]);
                            node3=node_overlying(node3,temp_value);
                        }
                        tree->insert(tree,f->id,f);
                        goto push;
                    }
                }
            }
        }
        vec.pop(&vec);
        // node=node_remove(node,node);
        // Node *node1=node;
        // node=(Node*)node->Next;
        // free(node1);
        //node1=NULL;
    }

    vec.clear(&vec);
    int_rb_tree_free(tree);
    // for(Node* nit=node3;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     hf=(template_hf*)(nit->value);
    //     f=hf->face;
    //     printf(" %d %d\n",f->vertices[0]->id,f->vertices[1]->id);
    // }
    //printf("end here\n");
    return node3;
}

//

static template_hf Mesh_opposite_halfface(template_hf hf)
{
    template_f*f=(template_f*)hf.face;
    if(f->halffaces[0]->cell==hf.cell)
    {
        return *(f->halffaces[1]);
    }
    else
    {
        return *(f->halffaces[0]);

    }
}

//


static iterator_v Mesh_v_beginn(struct Mesh * own)
{
    iterator_v it;
    iterator_v_init(&it);
    it.it=own->vertices->begin(own->vertices);
    return it;
}
static iterator_v Mesh_v_endn(struct Mesh*own)
{
   iterator_v it;
   iterator_v_init(&it);
   return it;
}

static iterator_hf Mesh_hf_beginn(struct  Mesh*own)
{
    iterator_hf it;
    iterator_hf_init(&it);
    it.it=own->halffaces->begin(own->halffaces);
    return it;
}
static iterator_hf Mesh_hf_endn(struct  Mesh*own)
{
    iterator_hf it;
    iterator_hf_init(&it);
    return it;
}

static iterator_f Mesh_f_beginn(struct  Mesh*own)
{
    iterator_f it;
    iterator_f_init(&it);
    it.it=own->faces->begin(own->faces);
    return it;
}
static iterator_f Mesh_f_endn(struct  Mesh*own)
{
    iterator_f it;
    iterator_f_init(&it);
    return it;
}

static iterator_c Mesh_c_beginn(struct Mesh*own)
{
    iterator_c it;
    iterator_c_init(&it);
    it.it=own->cells->begin(own->cells);
    return it;
}
static iterator_c Mesh_c_endn(struct Mesh*own)
{
     iterator_c it;
    iterator_c_init(&it);
    return it;
}

static iterator_v Mesh_fv_begin(struct Mesh* own,template_f f_)
{

    iterator_v it;
    iterator_v_init(&it);
    it.value=f_.vertices;
    return it;

}

static iterator_v Mesh_fv_end(struct Mesh* own,template_f f_)
{
    iterator_v it;
    iterator_v_init(&it);
    it.value=f_.vertices;
    it.i=f_.vertices_size;
    return it;

}
static iterator_v Mesh_hfv_begin(struct Mesh* own,template_hf hf)
{
    iterator_v it;
    iterator_v_init(&it);
    it.value=hf.vertices;
    return it;

}
static iterator_v Mesh_hfv_end(struct Mesh* own,template_hf hf)
{
    iterator_v it;
    iterator_v_init(&it);
    it.value=hf.vertices;
    it.i=hf.vertices_size;

    return it;
}
static iterator_v Mesh_cv_begin(struct Mesh* own,template_c c)
{
    iterator_v it;
    iterator_v_init(&it);
    it.value=c.vertices;
    return it;
}
static iterator_v Mesh_cv_end(struct Mesh* own,template_c c)
{
    iterator_v it;
    iterator_v_init(&it);
    it.value=c.vertices;
    it.i=c.vertices_size;
    return it;
}


static iterator_f Mesh_vf_begin(struct Mesh* own,template_v v)
{
    iterator_f iff;
    iterator_f_init(&iff);
    iff.vec=(v.facesn);

    // if(v.faces==NULL)
    // {
    //     return iff;
    // }
    // //***********
    // //可以通过node->traits记录起点增加速度
    // //***************
    // iff.node=*(v.faces);
    //iff.node=*(node_reverse(v.faces));

    return iff;

}
static iterator_f Mesh_vf_end(struct Mesh* own,template_v v)
{
    iterator_f iff;
    iterator_f_init(&iff);
    iff.vec=(v.facesn);
    iff.i=v.facesn.size;
    return iff;

}
iterator_hf Mesh_chf_begin(struct Mesh* own, template_c c)
{

    iterator_hf iff;
    iterator_hf_init(&iff);
    iff.value=c.halffaces;

    // if(c.halffaces==NULL)
    // {
    //     return iff;
    // }
    // // if(c.halffaces->traits==NULL)
    // // {
    // //     c.halffaces->traits=(void*)(node_reverse(c.halffaces));
    // // }
    // iff.node=*((Node*)(c.halffaces));

    return iff;
}
iterator_hf Mesh_chf_end(struct Mesh* own,template_c c)
{
    iterator_hf iff;
    iterator_hf_init(&iff);
    iff.value=c.halffaces;
    iff.i=c.halffaces_size;
    return iff;

}
static iterator_c Mesh_vc_begin(struct Mesh* own,template_v v)
{
    iterator_c it;
    iterator_c_init(&it);
    it.vec=(v.cellsn);
    // if(v.cells==NULL)
    // {
    //     return it;
    // }
    // it.node=*(v.cells);
    return it;

}
static iterator_c Mesh_vc_end(struct Mesh* own,template_v v)
{

    iterator_c it;
    iterator_c_init(&it);
    it.vec=(v.cellsn);
    it.i=v.cellsn.size;
    return it;
}
static iterator_edge Mesh_edge_begin(struct Mesh*own,template_v v,template_f f  )
{

    iterator_edge ie;
    iterator_edge_init(&ie);

    // for(Node* nit=v.faces;nit!=NULL;nit=(Node*)(nit->Next))
    // for(int i=0;i<v.facesn.size;i++)
    for(auto vfit=own->vf_begin(own,v);vfit!=own->vf_end(own,v);vfit++)
    {
        // template_f* f1=(template_f*)(nit->value);
        // template_f*f1=(template_f*)(v.facesn.extra_data[i]);
        template_f*f1=quote(vfit);
        if(f1->id==f.id)
        {
            ie.f=f1;
            break;
        }
    }
    for(int i=0;i<f.vertices_size;i++)
    {
        if(f.vertices[i]->id==v.id)
        {
            ie.v=f.vertices[i];
            break;
        }
    }
    ie.mark=ie.f;

    return ie;
}
Node* Mesh_vv_begin(struct Mesh* own,template_v v)
{
    // if(own->simplex!=1)
    // {
    //     printf("is not simplex\n");
    //     return NULL;
    // }
    Node* node=NULL;
    // printf("vv begin %d %d\n",v.id,v.facesn.size);
    iterator_f  myfit=Mesh_vf_begin(own,v);
    // printf("myfit.i :%d\n",myfit.i);
    iterator_f  endfit=Mesh_vf_end(own,v);
    // printf("endfit.i :%d\n",endfit.i);
    for(int i=0;i<v.facesn.size;i++)
    {
        template_f* f=(template_f*)(v.facesn.extra_data[i]);
        // printf("f id:%d\n",f->id);

    }

    for(iterator_f fit=Mesh_vf_begin(own,v);fit!=Mesh_vf_end(own,v);fit++)
    {
        // printf("fit.i:%d\n",fit.i);
        // if(fit.vec!=NULL)
        // {
        //     printf("shi \n");
        //     template_f* f=(template_f*)(v.facesn.extra_data[fit.i]);

        //     // template_f* f=(template_f*)(fit.vec->extra_data[fit.i]);
        //     printf("f id:%d\n",f->id);
        // }
        // f_.vec!=NULL?((template_f*)(f_.vec->extra_data[f_.i] ))
        // printf("vf id:%d\n",quote(fit)->id);
        for(auto vit=Mesh_fv_begin(own,*fit);vit!=Mesh_fv_end(own,*fit);vit++)
        {
            if((*vit).id!=v.id)
            {
                if(node_find(node,(void*)quote(vit))==NULL)
                {
                    node=node_overlying(node,(void*)quote(vit));
                }
            }
        }
    }
    return node;
}

LB_Vector_voidptr Mesh_vv_beginn(struct Mesh* own,template_v v)
{
    LB_Vector_voidptr vec;
    lb_vector_init_voidptr(&vec);
    // Node* node=NULL;
    for(iterator_f fit=Mesh_vf_begin(own,v);fit!=Mesh_vf_end(own,v);fit++)
    {
        for(auto vit=Mesh_fv_begin(own,*fit);vit!=Mesh_fv_end(own,*fit);vit++)
        {
            if((*vit).id!=v.id)
            {
                int id=-1;
                for(int i=0;i<vec.size;i++)
                {
                    if(*(vec.at(&vec,i))==quote(vit))
                    {
                        id=i;
                        break;
                    }
                }
                if(id==-1)
                {
                    vec.push_back(&vec,quote(vit));
                }
                // if(node_find(node,(void*)quote(vit))==NULL)
                // {
                //     // node=node_overlying(node,(void*)quote(vit));
                //     vec.push_back(&vec,quote(vit));
                // }
            }
        }
    }
    return vec;
}
//
Node* Mesh_vv_end(struct Mesh* own,template_v  vv)
{

    return NULL;
}

int Mesh_two_cell_is_connected(struct Mesh* own,template_c* c1,template_c* c2)
{
    bool re=false;
    if(c1==NULL||c2==NULL)
    {
        return re;
    }

    for(auto hfit=Mesh_chf_begin(own,*c1);hfit!=Mesh_chf_end(own,*c1);hfit++)
    {
        if(Mesh_s_opposite_halfface(quote(hfit))->cell==c2)
        {
            return true;
        }
    }
    return re;

}
static int Mesh_is_manifold_vertices(struct Mesh*own,template_v*v)
{
    // if(v->cells==NULL)
    if(v->cellsn.size<=0)
    {
        return  (own->dimension==0);
    }

    for(auto vfit=Mesh_vf_begin(own,*v);vfit!=Mesh_vf_end(own,*v);vfit++)
    {
        if(quote(vfit)->halffaces[0]->cell==NULL&&quote(vfit)->halffaces[1]->cell==NULL)
        {
            return 0;
        }

    }

    LB_Vector_voidptr vec;
    lb_vector_init_voidptr(&vec);
    vec.push_back(&vec,v->cellsn.extra_data[0]);

    LB_Vector1_voidptr vec1;
    lb_vector1_init_voidptr(&vec1);
    lb_vector1_resize_voidptr(&vec1,v->cellsn.size-1);
    memmove(vec1.extra_data,&(v->cellsn.extra_data[1]),sizeof(void*)*(v->cellsn.size-1) );
    int flag=0;
    do{
        flag=0;
        for(int i=0;i<vec1.size;i++)
        {
            template_c* c=(template_c*)(vec1.extra_data[i]);
            // for(Node* nnit=n1;nnit!=NULL;nnit=(Node*)(nnit->Next))
            for(int j=0;j<vec.size;j++)
            {
                if(Mesh_two_cell_is_connected(own,c,(template_c*)(*(vec.at(&vec,j)))))
                {
                    flag=1;
                    vec.push_back(&vec,c);
                    vec1.size--;
                    vec1.extra_data[i]=vec1.extra_data[vec1.size];
                    // lb_vector1_remove_voidptr(&vec1,c);
                    // n1=node_overlying(n1,nit->value);
                    // n2=node_remove(n2,nit);
                    // n2=node_delete_value(n2,nit->value);
                    break;
                }
            }
            if(flag==1)
            {
                break;
            }
        }

    }while(flag);
    int re=1;
    vec1.size==0?(re=1):(re=0);

    vec.clear(&vec);
    lb_vector1_clear_voidptr(&vec1);
    return re;

}
void Mesh_reset_c_prop(struct Mesh* mesh,int fr)
{
    for(auto cit=Mesh_c_beginn(mesh);cit!=Mesh_c_endn(mesh);cit++)
    {
        if(fr)
        {

            SAFE_FREE(quote(cit)->prop);
        }
        else
        {
            ((template_c*)(quote(cit)))->prop=NULL;
        }
    }
}
void Mesh_reset_v_prop(struct Mesh* mesh,int fr)
{
    // for(auto vit=mesh->vertices.begin();vit!=mesh->vertices.end();vit++)
    // {
    //     vit->second->prop=NULL;
    // }
    for(auto vit=Mesh_v_beginn(mesh);vit!=Mesh_v_endn(mesh);vit++)
    {
        if(fr)
        {
            SAFE_FREE(quote(vit)->prop);
        }
        else
        {
            ((template_v*)(quote(vit)))->prop=NULL;
        }
    }

}
void Mesh_reset_hf_prop(struct Mesh* mesh,int fr)
{
    for(auto hfit=Mesh_hf_beginn(mesh);hfit!=Mesh_hf_endn(mesh);hfit++)
    {
        if(fr)
        {
            SAFE_FREE(quote(hfit)->prop);
        }
        else
        {
            ((template_hf*)(quote(hfit)))->prop=NULL;
        }

    }
}
void Mesh_reset_f_prop(struct Mesh* mesh,int fr)
{
    for(auto fit=Mesh_f_beginn(mesh);fit!=Mesh_f_endn(mesh);fit++)
    {
        if(fr)
        {
            SAFE_FREE(quote(fit)->prop);
        }
        else
        {
            ((quote(fit)))->prop=NULL;
        }
    }
}

static void Mesh_printself(struct Mesh* own)
{
    printf("****************************\n\n");
    for(auto it=Mesh_v_beginn(own);it!=Mesh_v_endn(own);it++)
    {
        printf("id :%d \n",((template_v*)(quote(it)))->id);
        for(int j=0;j<((template_v*)(quote(it)))->point_size;j++)
        {
            printf("%lf  ",((template_v*)(quote(it)))->point[j]);
        }
        printf("\n");
    }

    printf("faces\n");
    for(auto it=Mesh_f_beginn(own);it!=Mesh_f_endn(own);it++)
    {
       printf("id:%d\n",((template_f*)(quote(it)))->id);
       if(((template_f*)(quote(it)))->prop!=NULL)
       {
            printf("prop ");
       }
       for(int j=0;j<((template_f*)(quote(it)))->vertices_size;j++)
       {
           printf("id :%d  ",((template_v*)(((template_f*)(quote(it)))->vertices[j]))->id);

       }
       printf("\n");
    }
    printf("cell\n");
    for(auto it=Mesh_c_beginn(own);it!=Mesh_c_endn(own);it++)
    {
        printf("cell id:%d\n",((template_c*)(quote(it)))->id);
        for(int j=0;j<((template_c*)(quote(it)))->vertices_size;j++)
        {
            printf("id %d  ",((template_v*)(((template_c*)(quote(it)))->vertices[j]))->id);
        }
        printf("\n");
    }
    // for(auto it=own->vertices.begin();it!=own->vertices.end();it++)
    // {
    //     printf("id :%d \n",it->second->id);
    //     for(int j=0;j<it->second->point_size;j++)
    //     {
    //         printf("%lf  ",it->second->point[j]);
    //     }
    //     printf("\n");
    // }
    // printf("faces\n");
    // for(auto it=own->faces.begin();it!=own->faces.end();it++)
    // {
    //    printf("id:%d\n",it->second->id);
    //    if(it->second->prop!=NULL)
    //    {
    //         printf("prop ");
    //    }
    //    for(int j=0;j<it->second->vertices_size;j++)
    //    {
    //        printf("id :%d  ",((template_v*)(it->second->vertices[j]))->id);

    //    }
    //    printf("\n");
    // }
    // printf("cell\n");
    // for(auto it=own->cells.begin();it!=own->cells.end();it++)
    // {
    //     printf("cell id:%d\n",it->second->id);
    //     for(int j=0;j<it->second->vertices_size;j++)
    //     {
    //         printf("id %d  ",((template_v*)(it->second->vertices[j]))->id);
    //     }
    //     printf("\n");
    // }
 printf("****************************\n\n");

}
Node* Mesh_intersection_two_faces(struct Mesh*own,template_f*f0,template_f*f1)
{
    Node * node=NULL;
    int flag=1;
    for(int i=0;i<f0->vertices_size;i++)
    {
        flag=1;
        for(int j=0;j<f1->vertices_size;j++)
        {
            if((f0->vertices[i])==f1->vertices[j])
            {
                flag=0;
                break;
            }
        }
        if(flag==0)
        {
            node=node_overlying(node,f0->vertices[i]);
        }
    }
    return node;
}


void Mesh_free(struct Mesh* own)
{
    if(own==NULL)
    {
        return;
    }
    SAFE_FREE(own->name);
    for(Node* nit=own->external_cells;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c*c=(template_c*)(nit->value);
        free_Cell(c);
    }
    free_node(own->external_cells);own->external_cells=NULL;
    // Mesh_external_cells_free_(own);
    for(auto iter=Mesh_c_beginn(own);iter!=Mesh_c_endn(own);iter++)
    {
        free_Cell((template_c*)(quote(iter)));
    }

    for(auto iter=Mesh_f_beginn(own);iter!=Mesh_f_endn(own);iter++)
    {
        free_Face((template_f*)(quote(iter)));
    }
    for(auto iter=Mesh_v_beginn(own);iter!=Mesh_v_endn(own);iter++)
    {
        free_Vertex((template_v*)(quote(iter)));
    }

    int_rb_tree_free(own->cells);
    int_rb_tree_free(own->faces);
    int_rb_tree_free(own->halffaces);
    int_rb_tree_free(own->vertices);
    own->cells=NULL;own->faces=NULL;own->halffaces=NULL;own->vertices=NULL;
    own->cell_id=0;own->vertex_id=0;own->face_id=0;own->halfface_id=0;


    // for(auto iter=own->cells.begin();iter!=own->cells.end();iter++)
    // {
    //     free_Cell(iter->second);
    // }
    // for(auto iter=own->faces.begin();iter!=own->faces.end();iter++)
    // {
    //     free_Face(iter->second);
    // }
    // for(auto iter=own->vertices.begin();iter!=own->vertices.end();iter++)
    // {
    //     free_Vertex(iter->second);
    // }
    // own->cells.clear();
    // own->faces.clear();
    // own->halffaces.clear();
    // own->vertices.clear();
}
//需要修改的函数

static Node* Mesh_non_manifold_vertices(struct Mesh* own)
{
    Node* re=NULL;
    for(auto iter=Mesh_v_beginn(own);iter!=Mesh_v_endn(own);iter++)
    {
        if(!Mesh_is_manifold_vertices(own,quote(iter)))
        {
            re=node_overlying(re,quote(iter));
        }

    }
    return re;
}

static Node* Mesh_one_dim_boundary(struct Mesh* own)
{
    Node* node=NULL;
    // for(auto it=own->faces.begin();it!=own->faces.end();it++)
    // {
    //     if(it->second->halffaces[0]->cell==NULL)
    //     {
    //         node=node_overlying(node,it->second->halffaces[0]);
    //     }
    //     else if(it->second->halffaces[1]->cell==NULL)
    //     {
    //         node=node_overlying(node,it->second->halffaces[1]);
    //     }
    // }
    for(auto it=Mesh_f_beginn(own);it!=Mesh_f_endn(own);it++)
    {
        if(((template_f*)(quote(it)))->halffaces[0]->cell==NULL)
        {
            node=node_overlying(node,((template_f*)(quote(it)))->halffaces[0]);
        }
        else if(((template_f*)(quote(it)))->halffaces[1]->cell==NULL)
        {
            node=node_overlying(node,((template_f*)(quote(it)))->halffaces[1]);
        }
    }


    return node;
}

void mesh_rearrange_elements_id(struct Mesh*own)
{
    int i=0;
    for(auto iter=Mesh_c_beginn(own);iter!=Mesh_c_endn(own);iter++)
    {
        quote(iter)->id=i;
        iter.it.it->key=i;
        i++;
    }
    i=0;
    for(auto iter=Mesh_f_beginn(own);iter!=Mesh_f_endn(own);iter++)
    {
        quote(iter)->id=i;
        iter.it.it->key=i;
        i++;
    }
    i=0;
    for(auto iter=Mesh_hf_beginn(own);iter!=Mesh_hf_endn(own);iter++)
    {
        quote(iter)->id=i;
        iter.it.it->key=i;
        i++;
    }
    i=0;
    for(auto iter=Mesh_v_beginn(own);iter!=Mesh_v_endn(own);iter++)
    {
        quote(iter)->id=i;
        iter.it.it->key=i;
        i++;
    }

}

// static void Mesh_external_cell_init_(struct Mesh* own)
// {
//     free_node(own->external_cell.halffaces);
//     own->external_cell.halffaces=NULL;
// //#ifdef MANIFOLD_REQUIRE
//     if(own->dimension==1&&own->manifold_require==1)
//     {
//           own->external_cell.halffaces= Mesh_one_dim_boundary(own);
// //            printf("external_cell size:%d\n",node_size(own->external_cell.halffaces));

//           return;
//     }
// //#endif
//     for(auto it=Mesh_f_beginn(own);it!=Mesh_f_endn(own);it++)
//     {
//         if(Mesh_nface_is_boundary(own,*((template_f*)(quote(it)))))
//         {
//             own->external_cell.halffaces=Mesh_node_of_boundary_face(own,((template_f*)(quote(it))));

//             break;
//         }
//     }
//     // for(auto it=own->faces.begin();it!=own->faces.end();it++)
//     // {
//     //     if(Mesh_nface_is_boundary(own,*(it->second)))
//     //     {
//     //         own->external_cell.halffaces=Mesh_node_of_boundary_face(own,it->second);
//     //         break;
//     //     }
//     // }
//     template_hf* hf1,*hf2;
//     Node* node=own->external_cell.halffaces;

//   //  printf("external_cell size:%d\n",node_size(node));
//     if(node==NULL)
//     {
//         return;
//     }

//     if(own->dimension==2)
//     {
//         own->external_cell.vertices_size=node_size(node);
//         SAFE_FREE(own->external_cell.vertices);
//         own->external_cell.vertices=(template_v**)malloc(sizeof(template_v*)*own->external_cell.vertices_size);

//     }

//     int i=0;
//     while(node!=NULL)
//     {
//         hf1=(template_hf*)(node->value);
//         if(hf1->vertices!=NULL)
//         {
//             free(hf1->vertices);
//             hf1->vertices=NULL;
//         }

//         hf2=Mesh_s_opposite_halfface(hf1);
//         hf1->vertices=(template_v**)malloc(sizeof(template_v*)*hf2->vertices_size);

//         hf1->vertices_size=hf2->vertices_size;
//         for(int i=0;i<hf2->vertices_size;i++)
//         {
//             hf1->vertices[i]=hf2->vertices[i];
// //              printf("external cell h_f id:%d ",((template_v*)(hf1->vertices[i]))->id);
//         }
//         if(hf1->vertices_size>=2)
//         {
//               template_v * temp=hf1->vertices[0];
//               hf1->vertices[0]=hf1->vertices[1];
//               hf1->vertices[1]=temp;
//          }
//          if(own->dimension==2)
//          {
//             own->external_cell.vertices[i]=hf1->vertices[0];
//          }
//          i++;
//          node=(Node*)(node->Next);
//     }
// }


static void my_fill_extern_cell_from_one_face(Mesh*own,template_f*f ,Int_RB_Tree* tree,template_c* c )
{
    Node* node =Mesh_node_of_boundary_face(own,f);
    if(node==NULL){return ;}
    for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf* hf=(template_hf*)(nit->value);
        tree->insert(tree,hf->face->id,hf->face);
    }
   // Node* node=c->halffaces;
    c->halffaces_size=node_size(node);
    c->halffaces=(template_hf**)malloc(sizeof(template_hf*)*c->halffaces_size);

    if(own->dimension==2)
    {
        c->vertices_size=node_size(node);
        //SAFE_FREE(own->external_cell.vertices);
        c->vertices=(template_v**)malloc(sizeof(template_v*)*c->vertices_size);
    }
    // printf("once external***********************\n");
    int i=0;
    for(Node* nit=node;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf*hf1=(template_hf*)(nit->value);
        c->halffaces[i]=hf1;
        if(hf1->vertices!=NULL)
        {
            free(hf1->vertices);
            hf1->vertices=NULL;
        }
        template_hf*hf2=Mesh_s_opposite_halfface(hf1);
        hf1->vertices=(template_v**)malloc(sizeof(template_v*)*hf2->vertices_size);
        hf1->vertices_size=hf2->vertices_size;
        memmove(hf1->vertices,hf2->vertices,sizeof(template_v*)*hf2->vertices_size);
        if(hf1->vertices_size>=2)
        {
              template_v * temp=hf1->vertices[0];
              hf1->vertices[0]=hf1->vertices[1];
              hf1->vertices[1]=temp;
        }
        // printf("external hf v id:%d %d  ** ",hf1->vertices[0]->id,hf1->vertices[1]->id);
        if(own->dimension==2)
        {
            c->vertices[i]=hf1->vertices[0];
            // printf("c vertices:%d ** \n",c->vertices[i]->id);
        }

        i++;
    }
}

static Node* mesh_get_all_external_cells(struct Mesh* own)
{
    Node* re=NULL;
    if(own->dimension==1)
    {
        //own->external_cell.halffaces= Mesh_one_dim_boundary(own);
        template_c* c=(template_c*)malloc(sizeof(template_c));
        Cell_init_(c);

        Node* halffaces=Mesh_one_dim_boundary(own);
        c->halffaces_size=node_size(halffaces);
        c->halffaces=(template_hf**)malloc(sizeof(template_hf*)*c->halffaces_size);
        int i=0;
        for(Node* nit=halffaces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            c->halffaces[i]=(template_hf*)(nit->value);
            i++;
        }

        free_node(halffaces);
        re=node_overlying(re,c);

//            printf("external_cell size:%d\n",node_size(own->external_cell.halffaces));
        return re;
    }
    Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree);

    for(auto it=Mesh_f_beginn(own);it!=Mesh_f_endn(own);it++)
    {
        if(Mesh_nface_is_boundary(own,*it)&&tree->find(tree,quote(it)->id )==NULL)
        {
            template_c* c=(template_c*)malloc(sizeof(template_c));
            Cell_init_(c);
            my_fill_extern_cell_from_one_face(own,quote(it),tree,c);
            re=node_overlying(re,c);

        }
    }

    int_rb_tree_free(tree);
    return re;
}


void Mesh_external_cells_init_(struct Mesh* own)
{
    for(Node* nit=own->external_cells;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c*c=(template_c*)(nit->value);
        free_Cell(c);
    }
    free_node(own->external_cells);
    own->external_cells=NULL;
    own->external_cells=mesh_get_all_external_cells(own);

}
void Mesh_external_cells_free_(struct Mesh*own)
{
    for(Node* nit=own->external_cells;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c*c=(template_c*)(nit->value);
        free_Cell(c);
    }
    free_node(own->external_cells);own->external_cells=NULL;

}
void Mesh_init(struct Mesh* own)
{
    own->name=NULL;
    own->faces=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    own->cells=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    own->halffaces=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    own->vertices=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init((own->cells));
    int_rb_tree_init((own->faces));
    int_rb_tree_init((own->vertices));
    int_rb_tree_init((own->halffaces));
    //own->cells.insert(&(own->cells),9,NULL);
    // own->vertices=(RB_Tree*)malloc(sizeof(RB_Tree));
    // rb_tree_init_int(own->vertices);
    // own->halffaces=(RB_Tree*)malloc(sizeof(RB_Tree));
    // rb_tree_init_int(own->halffaces);
    // own->faces=(RB_Tree*)malloc(sizeof(RB_Tree));
    // rb_tree_init_int(own->faces);
    // own->cells=(RB_Tree*)malloc(sizeof(RB_Tree));
    // rb_tree_init_int(own->cells);

    own->external_cells=NULL;
    //Cell_init_(&(own->external_cell));
    own->cell_id=0;own->vertex_id=0;own->face_id=0;own->halfface_id=0;
    own->get_vertexp=Mesh_get_vertexp;


    own->get_facep=Mesh_get_facep;
    own->get_halffacep=Mesh_get_halffacep;
    own->get_cellp=Mesh_get_cellp;
    own->get_vertex=Mesh_get_vertex;
    own->get_cell=Mesh_get_cell;
    own->get_face=Mesh_get_face;
    own->get_halfface=Mesh_get_halfface;
    own->get_facev=Mesh_get_facev;
    own->get_cellv=Mesh_get_cellv;
    own->get_cellf=Mesh_get_cellf;
    own->num_v=Mesh_num_v;
    own->num_c=Mesh_num_c;
    own->num_f=Mesh_num_f;
    own->num_hf=Mesh_num_hf;
    own->create_vertex=Mesh_create_vertex;
    own->create_vertexv=Mesh_create_vertexv;
    own->create_cell=Mesh_create_cell;
    own->create_face=Mesh_create_face;
    own->create_cellv=Mesh_create_cellv;
    own->create_cellf=Mesh_create_cellf;
    own->create_facev=Mesh_create_facev;
    own->create_halfface=Mesh_create_halfface;
    own->delete_cell=Mesh_delete_cell;
    own->delete_vertex=Mesh_delete_vertex;
    own->delete_face=Mesh_delete_face;
    own->delete_halfface=Mesh_delete_halfface;
    own->vertex_is_boundary=Mesh_vertex_is_boundary;
    own->face_is_boundary=Mesh_nface_is_boundary;
    own->cell_is_boundary=Mesh_cell_is_boundary;
    own->face_is_in_cell=Mesh_face_is_in_cell;
    own->vertex_is_in_cell=Mesh_vertex_is_in_cell;
    own->vertex_is_in_face=Mesh_vertex_is_in_face;

    own->node_of_boundary_face=Mesh_node_of_boundary_face;
    own->external_cells_init_=Mesh_external_cells_init_;
    own->external_cells_free_=Mesh_external_cells_free_;
    own->opposite_halfface=Mesh_opposite_halfface;
    own->s_opposite_halfface=Mesh_s_opposite_halfface;


    own->v_begin=Mesh_v_beginn;
    own->v_end=Mesh_v_endn;
    own->f_begin=Mesh_f_beginn;
    own->f_end=Mesh_f_endn;
    own->hf_begin=Mesh_hf_beginn;
    own->hf_end=Mesh_hf_endn;
    own->c_begin=Mesh_c_beginn;
    own->c_end=Mesh_c_endn;

    own->cv_begin=Mesh_cv_begin;
    own->cv_end=Mesh_cv_end;
    own->fv_begin=Mesh_fv_begin;
    own->fv_end=Mesh_fv_end;
    own->hfv_begin=Mesh_hfv_begin;
    own->hfv_end=Mesh_hfv_end;
    own->vc_begin=Mesh_vc_begin;
    own->vc_end=Mesh_vc_end;
    own->vf_begin=Mesh_vf_begin;
    own->vf_end=Mesh_vf_end;
    own->chf_begin=Mesh_chf_begin;
    own->chf_end=Mesh_chf_end;
    own->vv_begin=Mesh_vv_begin;
    own->vv_end=Mesh_vv_end;
    own->edge_begin=Mesh_edge_begin;
    own->intersection_two_faces=Mesh_intersection_two_faces;
    own->non_manifold_vertices=Mesh_non_manifold_vertices;
    own->is_manifold_vertices=Mesh_is_manifold_vertices;
    own->printself=Mesh_printself;

    own->reset_f_prop=Mesh_reset_f_prop;
    own->reset_hf_prop=Mesh_reset_hf_prop;
    own->reset_c_prop=Mesh_reset_c_prop;
    own->reset_v_prop=Mesh_reset_v_prop;


    own->simplex=1;
    own->dimension=-1;
    own->manifold_require=1;

    own->prop=NULL;
    meshT_init(&(own->traits));
}


Mesh* mesh_duplicate(struct Mesh*m  )
{
    Mesh* re=(Mesh*)malloc(sizeof(Mesh));
    Mesh_init(re);

    re->simplex=m->simplex;
    re->dimension=m->dimension;
    re->manifold_require=m->manifold_require;

    for(auto vit=m->v_begin(m);vit!=m->v_end(m);vit++ )
    {
        template_v* v=(template_v*)malloc(sizeof(template_v));
        Vertex_init_(v);
        v->id=quote(vit)->id;
        v->point_size=quote(vit)->point_size;
        v->point=(double*)malloc(sizeof(double)*v->point_size);
        memmove(v->point,quote(vit)->point,sizeof(double)*v->point_size);
        re->vertices->insert(re->vertices,v->id,v);
    }


    for(auto fit=m->f_begin(m);fit!=m->f_end(m);fit++)
    {
        template_v** vs=(template_v**)malloc(sizeof(template_v*)*quote(fit)->vertices_size);

        for(int i=0;i<quote(fit)->vertices_size;i++)
        {
            vs[i]=re->get_vertexp(re,quote(fit)->vertices[i]->id);
        }

        re->create_facev(re,vs,quote(fit)->vertices_size);
        free(vs);
    }

    for(auto cit=m->c_begin(m);cit!=m->c_end(m);cit++)
    {
        template_hf** hfs=(template_hf**)malloc(sizeof(template_hf*)*quote(cit)->halffaces_size);

        int i=0;
        for(auto chf=m->chf_begin(m,*cit);chf!=m->chf_end(m,*cit);chf++)
        {
            template_v** vs=(template_v**)malloc(sizeof(template_v*)*quote(chf)->vertices_size);

            for(int j=0;j<quote(chf)->vertices_size;j++)
            {
                vs[j]=re->get_vertexp(re,quote(chf)->vertices[j]->id) ;
            }
            template_f* f= re->get_facev(re,vs,quote(chf)->vertices_size);
            template_hf*hf=re->create_halfface(re,f,vs,quote(chf)->vertices_size);

            hfs[i]=hf;
            free(vs);
            i++;
        }

        re->create_cellf(re,hfs,quote(cit)->halffaces_size );
        free(hfs);
    }

    return re;
}
